## A few notes:
## 1. run_profile is my profile function so I can determine where
## 2. run_admb can run profile as well, but parameters difficult to control
## 3. perhaps most important, I turn off and on parameters, by removing the 
##    the "init_" out in the PARAMETER_SECTION and dumping them into the
##    PRELIMINARY_CALCS_SECTION; meaning, the parameters in the original .tpl all
##    need to be inside the PARAMETER_SECTION.
## 4. the format for the profile parameter is: "p_ORIGINALPARAMETERNAME".
## 5. variable declaration NEEDS TO BE IN ORDER
## My ADMB directory: /opt/local/admb/

write_admb_data <- function(filename, dat) {
  filename <- paste(filename, ".dat", sep = "")
  for (i in 1:length(dat)) {
    write(paste("#", names(dat)[i]), filename, append = ifelse(i == 1, F, T))
    if (is.null(dim(dat[[i]])))
      write(dat[[i]], filename, ncol = length(dat[[i]]), append = T)
    if (!is.null(dim(dat[[i]])))
      write(t(dat[[i]]), filename, ncol = ncol(dat[[i]]), append = T)
  }
}

write_admb_inits <- function(filename, inits) {
  filename <- paste(filename, ".pin", sep = "")
  for (i in 1:length(inits)) {
    write(paste("#", names(inits)[i]), filename, append = ifelse(i == 1, F, T))
    write(inits[[i]], filename, ncol = length(inits[[i]]), append = T)
  }
}

read_admb_results <- function(filename) {
  out <- list()
  std <- read_sd_results(filename)
  filename <- paste(filename, ".rep", sep = "")
  out.read <- readLines(filename)
  out.grep <- grep("# ", out.read)
  out.names <- gsub("# ", "", out.read[out.grep])
  for (i in 1:length(out.names)) {
    if (i < length(out.names))
      ind <- (out.grep[i]+1):(out.grep[i+1]-1)
    if (i == length(out.names))
      ind <- (out.grep[i]+1):length(out.read)
    out[[i]] <- matrix(na.omit(as.numeric(unlist(strsplit(out.read[ind], " ")))),
                       ncol = length(ind))
    if (length(ind) == 1)
      out[[i]] <- as.vector(out[[i]])
  }
  names(out) <- out.names
  out <- c(out, list(std = std))
  out
}

read_profile_results <- function(paramname) {
  out <- list() # NOTE: need "p_" in front of likeprof parameter
  filename <- paste("", paramname, ".plt", sep = "")
  out.read <- readLines(filename)
  out.grep <- grep("[A-Z]", out.read)
  out.len <- max(diff(out.grep)) - 2
  out.proflike <- grep("Profile likelihood", out.read) + 1
  out.normapprox <- grep("Normal approximation", out.read)[1] + 1
  out$proflike <- t(sapply(out.read[out.proflike:(out.proflike+out.len)], function(x)
                           as.numeric(strsplit(x, " ")[[1]][-1]), USE.NAMES = F))
  out$normapprox <- t(sapply(out.read[out.normapprox:(out.normapprox+out.len)], function(x)
                             as.numeric(strsplit(x, " ")[[1]][-1]), USE.NAMES = F))
  out
} 

read_sd_results <- function(filename) {
  out <- list()
  filename <- paste(filename, ".std", sep = "")
  if (file.exists(filename)) {
    if (file.info(filename)$size > 0) {
      out.read <- read.table(filename, skip = 1)
      names(out.read) <- c("index", "name", "value", "std dev")
      return(out.read)
    }
    else
      NULL
  }
  else
    NULL
}

restrict_admb_params <- function(filename, fixed) {
  file.tpl <- readLines(paste(filename, ".tpl", sep = ""))
  file.sec <- which(file.tpl == "PARAMETER_SECTION" |
                    file.tpl == "PRELIMINARY_CALCS_SECTION")
  file.fix <- file.tpl[file.sec[1]:file.sec[2]]
  for (i in 1:length(fixed)) {
    ind <- grep(paste("init_.*[[:space:]]", names(fixed)[i], sep = ""), file.fix)
    file.tpl[file.sec[1]-1+ind] <- gsub("init_", "", file.fix[ind])
    file.tpl <- append(file.tpl, paste(" ", names(fixed)[i], "=", fixed[i], ";"),
                       file.sec[2])
  }
  write(file.tpl, paste(filename, "_", paste(names(fixed), collapse = ""), ".tpl", sep = ""))
}

restrict_admb_params2 <- function(filename, fixed) {
  file.tpl <- readLines(paste(filename, ".tpl", sep = ""))
  file.sec <- which(file.tpl == "DATA_SECTION" |
                    file.tpl == "PARAMETER_SECTION")
  file.sec <- c(file.sec, grep("SECTION", file.tpl)
                [which(file.sec[2] == grep("SECTION", file.tpl)) + 1])
  file.fix <- file.tpl[file.sec[2]:file.sec[3]]
  for (i in 1:length(fixed)) {
    ind <- grep(paste("init_.*[[:space:]]", names(fixed)[i], sep = ""), file.fix)+i
    file.tpl[file.sec[2]-1+ind] <- gsub(file.fix[ind], "", file.fix[ind])
    file.tpl <- append(file.tpl, file.fix[ind], file.sec[2]-i)
  }
  write(file.tpl, paste(filename, "_", paste(names(fixed), collapse = ""), ".tpl", sep = ""))
}

run_admb <- function(filename, dat, init, fixed = NULL, re = FALSE,
                     is = F, is.N = 100,
                     profile = FALSE, prof_param, overwrite = FALSE,
                     clean = c("all","sys","output")) {
  clean_admb(filename, which = clean[1])
  if (length(fixed) > 0) {
    restrict_admb_params(filename, fixed)
    filename <- paste(filename, "_", paste(names(fixed), collapse = ""), sep = "")
  }
  write_admb_data(filename, dat)
  write_admb_inits(filename, init)
  if (overwrite == F) {
    out.compile <- system(paste("admb", ifelse(re, " -r", ""), " ",
                                filename, sep = ""),
                          intern = TRUE)
    if (length(grep("error", out.compile)) > 0)
      stop(out.compile)
  }
  out.run <- system(paste(ifelse(.Platform$OS.type == "unix", "./", "")
                          , filename, " ",
                          ifelse(re & is, "-is", ""), " ",
                          ifelse(re & is, is.N, ""), " ",                          
                          ifelse(profile, " -lprof", ""),
                          sep = ""), intern = TRUE)
  if (length(grep("Warning", out.run)) > 0)
    print(out.run[grep("Warning", out.run)])
  if (length(grep("Error", out.run)) > 0)
    print(out.run[grep("Error", out.run)])
  if (length(grep("out of bounds", out.run)) > 0)
    print(out.run[grep("out of bounds", out.run)])
  posdef = TRUE
  if (length(grep("Hessian does not appear to be positive definite", out.run)) > 0)
    posdef = FALSE
  if (length(grep("Can't make a positive definite Hessian", out.run)) > 0)
    posdef = FALSE
  if (length(grep("Error matrix not positive definite in choleski_decomp", out.run)) > 0)
    posdef = FALSE
  if (length(grep("Found positive definite Hessian with trust region method", out.run)) > 0)
    posdef = FALSE
  if (length(grep("Estimated covariance matrix may not be positive definite", out.run)) > 0)
    posdef = FALSE
  out <- read_admb_results(filename)
  if (file.exists(paste(filename, ".std", sep = "")))
    out <- c(out, list(std = read_sd_results(filename)))
  if (profile == T)
    out <- c(out, list(profile = read_profile_results(prof_param)))
  if (file.exists(paste(filename, ".cor", sep = "")))
    out <- c(out, read_cor_file(filename))
  c(out, posdef = posdef)
}

run_profile <- function(filename, dat, init, name, values, overwrite = F) {
  out <- rep(0, length(values))
  fixed <- values; names(fixed) <- rep(name, length(values))
  for (i in 1:length(values)) 
    out[i] <- run_admb(filename, dat, init, fixed[i], overwrite = overwrite)$obj_fun
  cbind(values = values, out = out)
}

# from Ben Bolker
clean_admb <- function(filename,which=c("all","sys","output"),profpars=NULL) {
  tplfile <- paste(filename,"tpl",sep=".")
  ## need logic here!  enumeration of all files ???
  sys.ext <- c("htp","cpp","o","rhes","bar","eva",
               "bgs","ecm","luu","mc2","mcm","tpl.bak")
  input.ext <- c("pin","dat")
  output.ext <- c("log","cor","std", "par","psv","hst","prf")
  more.ext <- c("b01","b02","obj","p01","p02","r01","r02")
  sys2.ext <- c("out","cout")
  other <- c("eigv.rpt","fmin.log","variance","sims",
             "hesscheck","hessian.bin","dgs2","diags",
             paste("admodel",c("dep","hes","cov"),sep="."))
  which <- match.arg(which)
  if (which=="all") {
    ## erase only targeted extensions -- NOT everything with the basename except the tpl file
    delfiles <-  paste(filename,c(sys.ext,input.ext,output.ext,sys2.ext,more.ext),sep=".")
    ## list.files(pattern=paste("^",filename,"\\..*",sep=""))
    ##    delfiles <- setdiff(delfiles,tplfile)
    delfiles <- c(delfiles)
    if (!is.null(profpars)) {
      delfiles <- c(delfiles,paste(profpars,".plt",sep=""))
    }
  } else {
    stop("only 'all' option is currently implemented")
  } 
  unlink(delfiles)
}

read_mcmc_results <- function(filename, inits) {
  file.names <- rep(names(inits), unlist(lapply(inits, length)))
  filen <- file(paste(filename, ".psv", sep = ""), "rb")
  nopar <- readBin(filen, what = integer(), n = 1)
  mcmc <- readBin(filen, what = numeric(), n = nopar * 10000)
  mcmc <- matrix(mcmc, byrow = TRUE, ncol = nopar)
  mcmc <- as.data.frame(mcmc)
  names(mcmc) <- file.names
  close(filen)
  mcmc
}

read_cor_file <- function(filename) {
  filename <- paste(filename, ".cor", sep = "")
  hess <- readLines(filename)
  file.len <- length(hess)
  hess <- strsplit(hess[1], " ")[[1]]
  hess <- as.numeric(hess[length(hess)])
  vals <- read.table(filename, skip = 2, fill = TRUE, head = FALSE, col.names = c("index", "name", "value", "std_dev", 1:(file.len-2)))
  list(hess = hess, cor = vals)
}
